#!/usr/bin/Rscript
setwd("/home/zhang/BMK/test/cuffdiff")

# 2016.11.16
# setwd("/home/zhang/Summer/cerna结果")
source("https://www.bioconductor.org/biocLite.R")
library(clusterProfiler)
library(openxlsx)
library(org.Hs.eg.db)
library(RDAVIDWebService)
library(DOSE)
library(igraph)
library(KEGG.db)

# 做do疾病富集注释，效果极差
do.do <- function(genelist){
  do <- enrichDO(genelist$ENTREZID, pvalueCutoff = 0.8, pAdjustMethod = "BH")
  do <- setReadable(do, 'org.Hs.eg.db')
  do <- summary(do)
  return(do)
}
# 做NCG的疾病富集注释，效果极差
do.ncg <- function(genelist){
  ncg <- enrichNCG(genelist$ENTREZID, pvalueCutoff = 0.8, pAdjustMethod = "BH")
  ncg <- setReadable(ncg, 'org.Hs.eg.db')
  ncg <- summary(ncg)
  return(ncg)
}
# 做kegg通路富集注释，kegg无法使用setReadable，所以ID没法转化为gene name这个很不爽
do.kegg <- function(genelist){
  kegg <- enrichKEGG(genelist$ENTREZID, pvalueCutoff = 0.05, pAdjustMethod = "BH")
  # kegg <- setReadable(kegg, 'org.Hs.eg.db')
  kegg <- summary(kegg)
  return(kegg)
}
# 做 GO Term的注释
do.go <- function(gene.df){
  real.do <- function(gene.df, type){
    cat(paste(as.character(type), sep = "\n"))
    ego2 <- enrichGO(gene         = gene.df$ENSEMBL,
                     OrgDb         = org.Hs.eg.db,
                     keytype       = 'ENSEMBL',
                     ont           = as.character(type),
                     pAdjustMethod = "BH",
                     pvalueCutoff  = 0.01,
                     qvalueCutoff  = 0.05)
    ego2 <- setReadable(ego2, OrgDb = org.Hs.eg.db)
    return(summary(ego2))
  }

  cc <- real.do(gene.df, "CC")
  bp <- real.do(gene.df, "BP")
  mf <- real.do(gene.df, "MF")

  if(nrow(cc) == 0){
    stop("cc没有富集")
  }else if(nrow(bp) == 0){
    stop("bp没有富集")
  }else if(nrow(mf) == 0){
    stop("mf没有富集")
  }
  
  cc$type <- "cellular component"
  bp$type <- "biological process"
  mf$type <- "molecular funtion"

  all <- rbind(cc, bp)
  all <- rbind(all, mf)

  return(all)
}

# 做结果的网络BC值
do.bc <- function(miranda){
  g <- graph.data.frame(miranda[,1:2], directed = T)

  # normalized = FALSE
  BC <- as.data.frame(betweenness(g, directed = F, normalized = F))
  colnames(BC)[1] <- "BC_not"

  # normalized = TRUE
  BC1 <- as.data.frame(betweenness(g, directed = F, normalized = T))
  colnames(BC1)[1] <- "BC_normalized"

  BC <- cbind(BC, BC1)

  BC$type <- "microRNA"

  mrna <- as.character(unique(miranda[miranda[,3] == "mrna",2]))
  BC[mrna,"type"] <- "mrna"

  lnc <- as.character(unique(miranda[miranda[,3] == "lnc",2]))
  BC[lnc,"type"] <- "lnc"

  BC$RNA <- rownames(BC)

  return(BC[,c(4,1:3)])
}

# 处理所有差异表达的mrna的结果
mrna <- read.table("mrna/gene_exp.diff", header = T)
mrna <- mrna[mrna$significant == "yes", ]
mrna <- as.character(unique(mrna$gene))
# 然后全套大保健
eg <- bitr(mrna, fromType = "SYMBOL", toType = c("ENTREZID", "ENSEMBL"), OrgDb = "org.Hs.eg.db")

do <- do.do(eg)
ncg <- do.ncg(eg)
kegg <- do.kegg(eg)
go <- do.go(eg)

wb = createWorkbook()
addWorksheet(wb, "KEGG")
addWorksheet(wb, "GO Tearm")
addWorksheet(wb, "DO")
addWorksheet(wb, "NCG")
addWorksheet(wb, "plog")
writeData(wb,1, kegg)
writeData(wb,2, go)
writeData(wb,3, do)
writeData(wb,4, ncg)
writeData(wb,5, eg)
saveWorkbook(wb, "mrna_david.xlsx")


# cerna中mRNA的通路注释
mrna <- read.table("mrna.cerna")

eg <- bitr(mrna[,1], fromType = "SYMBOL", toType = c("ENTREZID", "ENSEMBL"), OrgDb = "org.Hs.eg.db")

do <- do.do(eg)
ncg <- do.ncg(eg)
kegg <- do.kegg(eg)
go <- do.go(eg)


wb = createWorkbook()
addWorksheet(wb, "KEGG")
addWorksheet(wb, "GO Tearm")
addWorksheet(wb, "DO")
addWorksheet(wb, "NCG")
addWorksheet(wb, "plog")
writeData(wb,1, kegg)
writeData(wb,2, go)
writeData(wb,3, do)
writeData(wb,4, ncg)
writeData(wb,5, eg)
saveWorkbook(wb, "cenra_mrna_david.xlsx")

# 所有lnc上下游300kb靶基因的各种通路注释
mrna <- read.table("/home/zhang/BMK/test/cuffdiff/lnctarget/lnc.mrnalist")

eg <- bitr(mrna[,1], fromType = "SYMBOL", toType = c("ENTREZID", "ENSEMBL"), OrgDb = "org.Hs.eg.db")

do <- do.do(eg)
ncg <- do.ncg(eg)
kegg <- do.kegg(eg)
go <- do.go(eg)

wb = createWorkbook()
addWorksheet(wb, "KEGG")
addWorksheet(wb, "GO Tearm")
addWorksheet(wb, "DO")
addWorksheet(wb, "NCG")
addWorksheet(wb, "plog")
writeData(wb,1, kegg)
writeData(wb,2, go)
writeData(wb,3, do)
writeData(wb,4, ncg)
writeData(wb,5, eg)
saveWorkbook(wb, "lnc_david.xlsx")

# 处理cerna中lnc的靶基因的结果的
lnc <- read.table("/home/zhang/BMK/test/cuffdiff/lnctarget/lnc_cerna.mrnalist")
# 然后全套大保健
eg <- bitr(lnc[,1], fromType = "SYMBOL", toType = c("ENTREZID", "ENSEMBL"), OrgDb = "org.Hs.eg.db")

do <- do.do(eg)
ncg <- do.ncg(eg)
kegg <- do.kegg(eg)
go <- do.go(eg)


wb = createWorkbook()
addWorksheet(wb, "KEGG")
addWorksheet(wb, "GO Tearm")
addWorksheet(wb, "DO")
addWorksheet(wb, "NCG")
addWorksheet(wb, "plog")
writeData(wb,1, kegg)
writeData(wb,2, go)
writeData(wb,3, do)
writeData(wb,4, ncg)
writeData(wb,5, eg)
saveWorkbook(wb, "cenra_lnc_david.xlsx")


# 先做BC值的统计
# 两个输入文件就是miranda的输出结果，提交给cerna.jar的那个输入文件
# lnc和mrna分开，全部经过ago验证的那些
cerna <- read.csv("cerna结果/cerna_diff.csv")
cerna1 <- cerna[,c("RNA_1", "gene_name_1"), drop = F]
cerna2 <- cerna[,c("RNA_2", "gene_name_2"), drop = F]
colnames(cerna2) <- colnames(cerna1)
cerna <- unique(rbind(cerna1, cerna2))


mrna <- read.table("miranda结果/mrna.test")
lnc <- read.table("miranda结果/lnc.test")


mrna$V3 <- "mrna"
mrna <- merge(mrna, cerna, by.x = "V2", by.y = "RNA_1")
mrna <- mrna[,c(2,4,3)]

lnc$V3 <- "lnc"
lnc <- merge(lnc, cerna, by.x = "V2", by.y = "RNA_1")
lnc <- lnc[,c(2,4,3)]

miranda <- unique(rbind(mrna, lnc))

bc <- do.bc(miranda)
diff <- bc[bc$BC_normalized >= 0.02,]

cyto <- merge(miranda, diff, by.x = "V1", by.y = "RNA")
cyto <- cyto[,1:3]


wb = createWorkbook()
addWorksheet(wb, "miranda")
addWorksheet(wb, "BC")
addWorksheet(wb, "BCBigThan0.02")
addWorksheet(wb, "cyto")
writeData(wb,1, miranda)
writeData(wb,2, bc)
writeData(wb,3, diff)
writeData(wb,4, cyto)
saveWorkbook(wb, "mirandaBC_diff.xlsx")


#####  做superdiff

cerna <- read.xlsx2("cerna结果/cerna_superdiff.xlsx", 1)
cerna1 <- cerna[,c("RNA_1", "gene_name_1"), drop = F]
cerna2 <- cerna[,c("RNA_2", "gene_name_2"), drop = F]
colnames(cerna2) <- colnames(cerna1)
cerna <- unique(rbind(cerna1, cerna2))


mrna <- read.table("miranda结果/mrna.test")
lnc <- read.table("miranda结果/lnc.test")


mrna$V3 <- "mrna"
mrna <- merge(mrna, cerna, by.x = "V2", by.y = "RNA_1")
mrna <- mrna[,c(2,4,3)]

lnc$V3 <- "lnc"
lnc <- merge(lnc, cerna, by.x = "V2", by.y = "RNA_1")
lnc <- lnc[,c(2,4,3)]

miranda <- unique(rbind(mrna, lnc))

bc <- do.bc(miranda)
diff <- bc[bc$BC_normalized >= 0.02,]

cyto <- merge(miranda, diff, by.x = "V1", by.y = "RNA")
cyto <- cyto[,1:3]


wb = createWorkbook()
addWorksheet(wb, "miranda")
addWorksheet(wb, "BC")
addWorksheet(wb, "BCBigThan0.02")
addWorksheet(wb, "cyto")
writeData(wb,1, miranda)
writeData(wb,2, bc)
writeData(wb,3, diff)
writeData(wb,4, cyto)
saveWorkbook(wb, "mirandaBC_superdiff.xlsx")
